/*

This file generates Figure 3 in the main text.

*/


*>--------- define rounding function for regs---------------------<
* copied from rounding_regs.do Census file provided

capture program drop round
program define round, eclass
   mat A = e(b)
   forval j = 1/`=colsof(matrix(A))' { 
      mat A[1,`j'] = round(A[1,`j'],10^(floor(log10(abs(A[1,`j'])))-3))
   } 
   ereturn repost b = A

   if e(N)>=15&e(N)<=99 ereturn scalar N = round(e(N), 10)
   if e(N)>=100&e(N)<=999 ereturn scalar N = round(e(N), 50)
   if e(N)>=1000&e(N)<=9999 ereturn scalar N = round(e(N), 100)
   if e(N)>=10000&e(N)<=99999 ereturn scalar N = round(e(N), 500)
   if e(N)>=100000&e(N)<=999999 ereturn scalar N = round(e(N), 1000)
   if e(N)>=1000000 ereturn scalar N = round(e(N),10^(floor(log10(e(N)))-3))

   if e(N)<15 {
   scalar N_less15=1
   ereturn scalar N=.
   }

   if e(df_r)>=15&e(df_r)<=99 ereturn scalar df_r = round(e(df_r), 10)
   if e(df_r)>=100&e(df_r)<=999 ereturn scalar df_r = round(e(df_r), 50)
   if e(df_r)>=1000&e(df_r)<=9999 ereturn scalar df_r = round(e(df_r), 100)
   if e(df_r)>=10000&e(df_r)<=99999 ereturn scalar df_r = round(e(df_r), 500)
   if e(df_r)>=100000&e(df_r)<=999999 ereturn scalar df_r = round(e(df_r), 1000)
   if e(df_r)>=1000000 ereturn scalar df_r = round(e(df_r),10^(floor(log10(e(df_r)))-3))

   if e(df_r)<15 {
   scalar df_less15=1
   ereturn scalar df_r=.
   }


   capture ereturn scalar r2 = round(e(r2),10^(floor(log10(e(r2)))-3))
   capture ereturn scalar r2_a = round(e(r2_a),10^(floor(log10(e(r2_a)))-3))
   capture ereturn scalar F = round(e(F),10^(floor(log10(e(F)))-3))
   capture ereturn scalar rmse = round(e(rmse),10^(floor(log10(e(rmse)))-3))
   capture ereturn scalar mss = round(e(mss),10^(floor(log10(e(mss)))-3))
   capture ereturn scalar rss = round(e(rss),10^(floor(log10(e(rss)))-3))
   capture ereturn scalar rkf = round(e(rkf),10^(floor(log10(e(rkf)))-3))

end


*>---------------------------- electricity intensity --------------------------<

use "/projects/data/dataSTATA/combined/combined_productivity_energy_4digit_4m.dta", clear

ivreghdfe lelec_int  (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst)
keep if e(sample)

ivreghdfe lco2_int_i_both  (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst)
keep if e(sample)

ivreghdfe lbtu_int_i_both  (lelec_price lelec_price_init = instr_B_cl instr_B_ng instr_B_pa  instr_B_cl_init instr_B_ng_init instr_B_pa_init) [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst)
keep if e(sample)

reghdfe lelec_price_init instr_B_cl instr_B_ng instr_B_pa instr_B_cl_init instr_B_ng_init instr_B_pa_init   [pw=wt_cmf], absorb(i.y##i.firstyear##i.bestnaics i.fipsst##i.bestnaics) cluster(fipsst)
keep if e(sample)

* round variables before taking the average
do /projects/programs/codeSTATA/B_analysis/rounding_4sigdig_v2 firstyear year ee pe elec_price_init tvs co2_int_i_both btu_int_i_both
gen  elec_price_round = ee / pe 
gen  elec_int_round = pe / tvs


preserve

collapse (mean) elec_int_round co2_int_i_both btu_int_i_both oe oe_scale oh [pw = wt_cmf], by(year)
ren elec_int_round elec_int

replace btu = btu * 100
la var elec_int "Electricity"
la var co2_int_i_both "CO2"
la var btu_int_i_both "BTU"


* measure change relative to 1976
foreach var in elec_int co2_int_i_both btu_int_i_both oe oe_scale oh {
	gen `var'76 = `var' if year == 1976
	egen `var'_1976 = mean(`var'76)
	replace `var' = `var' / `var'_1976
}

* round before graphing
do /projects/programs/codeSTATA/B_analysis/rounding_4sigdig_v2 elec_int co2_int_i_both btu_int_i_both oe oe_scale

tw (line co2_int_i_both year, lpattern(dash) lcolor(cranberry)) (line elec_int year, lpattern(solid) lcolor(navy)) (line btu_int_i_both year, lpattern(shortdash) lcolor(green)), graphr(color(white) lwidth(vthick)) yscale(noline) xtit("Year") subtit("") ytitle("Energy Intensity") ylabel(0(0.5)2)
graph export "/projects/results/figures/intensity_trends.png", as(png) name("Graph") replace



*>--------------------------- electricity productivity  -----------------------<

la var oe_scale "Relative Energy Productivity"
la var oh "Total Factor Productivity"

tw (line oe_scale year, lpattern(solid) lcolor(navy)) (line oh year, lpattern(dash) lcolor(cranberry)), graphr(color(white) lwidth(vthick)) yscale(noline) xtit("Year") subtit("") ytitle("Productivity") ylabel(0.5(0.5)3.0) 
graph export "/projects/results/figures/productivity_trends.png", as(png) name("Graph") replace

restore

* export data to support file for disclosure review
export delimited using "/projects/disclosure/20210715/support/graphs_data.csv", replace
